<!DOCTYPE html>
<html>
<head>
	<meta charset="utf-8">
	<meta http-equiv="X-UA-Compatible" content="IE=edge">
	<meta name="viewport" content="width=device-width, initial-scale=1, user-scalable=no">
	<link rel="stylesheet" href="/assets/css/atom-one-light.css">
    
        <title>加速Python</title>
		<link rel="stylesheet" type="text/css" href="/assets/css/002.css">
    
	<link rel="stylesheet" href="/assets/css/font-awesome.min.css">
	<link rel="shortcut icon" href="/assets/img/favicon.ico" type="image/x-icon">
	<link rel="icon" href="/assets/img/favicon.ico" type="image/x-icon">
	<script src="/assets/js/highlight.pack.js"></script>
	<script>hljs.initHighlightingOnLoad();</script>

	

	<script async src="https://www.googletagmanager.com/gtag/js?id=UA-140127665-1"></script>
<script>
  window.dataLayer = window.dataLayer || [];
  function gtag(){dataLayer.push(arguments);}
  gtag('js', new Date());

  gtag('config', 'UA-140127665-1');
</script>


</head>
<body>
	<div class="wrapper">
		<div class="default_title">
			<img src="/assets/img/mycomputer.png" />
			
				<h1>NaNg's blog</h1>
			
		</div>
		<ul class="topbar">
	<a href="/pages/me.html"><li><u>A</u>bout</li></a>
	<a href="/pages/links.html"><li><u>L</u>inks</li></a>
	<a href="https://www.dropbox.com/sh/nhy3m3cvojizkk5/AABH8mt5gh3PiBrtWnCxE87ma?dl=0"><li><u>D</u>ropbox</li></a>
</ul>
		<div class="tag_list">
			<ul id="tag-list">
				<li><a href="/" ><img src="/assets/img/disk.png" />(C:)</a>
			<ul>
				
				
				<li><a href="/tag/3d/" title="3d"><img src="/assets/img/folder.ico" />3d</a></li>
				
				<li><a href="/tag/bioinformatics/" title="bioinformatics"><img src="/assets/img/folder.ico" />bioinformatics</a></li>
				
				<li><a href="/tag/notes/" title="notes"><img src="/assets/img/folder.ico" />notes</a></li>
				
				<li><a href="/tag/others/" title="others"><img src="/assets/img/folder.ico" />others</a></li>
				
				<li><a href="/tag/sci-fiction/" title="sci fiction"><img src="/assets/img/folder.ico" />sci fiction</a></li>
				
			</ul>
				</li>
			</ul>
		</div>
		<div class="post_list">
			
				<ul>
					
					<li>
						<a href="https://ani-net-project.gitee.io/index.html" title="AniNet">
								<img class="small-icon" src="/assets/img/aninet.png" title="AniNet" />AniNet
						</a>
					</li>
					
					<li>
						<a href="/examples/boids/index.html" title="Boids">
								<img class="small-icon" src="/assets/img/bird.png" title="Boids" />Boids
						</a>
					</li>
					
					<li>
						<a href="/examples/EM-alg/1_how_it_works.html" title="EM 算法 1. how it works">
								<img class="small-icon" src="/assets/img/notebook.ico" title="EM 算法 1. how it works" />EM 算法 1. how it works
						</a>
					</li>
					
					<li>
						<a href="/examples/pyvm/pyvm_ch0-3_cn.html" title="翻译：Inside Python Virtual Machine（前三章）">
								<img class="small-icon" src="/assets/img/html_ie.ico" title="翻译：Inside Python Virtual Machine（前三章）" />翻译：Inside Python Virtual Machine（前三章）
						</a>
					</li>
					
					<li>
						<a href="/20190624/cLife.html" title="器官工业幻想">
								<img class="small-icon" src="/assets/img/file.ico" title="器官工业幻想" />器官工业幻想
						</a>
					</li>
					
					<li>
						<a href="/examples/pubnet/Network_statistic.html" title="PubNet network statistic">
								<img class="small-icon" src="/assets/img/notebook.ico" title="PubNet network statistic" />PubNet network statistic
						</a>
					</li>
					
					<li>
						<a href="/examples/pubnet/sample.html" title="PubMed bio-conception network example">
								<img class="small-icon" src="/assets/img/net.png" title="PubMed bio-conception network example" />PubMed bio-conception network example
						</a>
					</li>
					
					<li>
						<a href="/examples/hpo_enrich/example_sagd_00055.html" title="HPO enrichment example">
								<img class="small-icon" src="/assets/img/notebook.ico" title="HPO enrichment example" />HPO enrichment example
						</a>
					</li>
					
					<li>
						<a href="/20190222/KS-Test.html" title="Two sample Kolmogorov-Smirnov test">
								<img class="small-icon" src="/assets/img/file.ico" title="Two sample Kolmogorov-Smirnov test" />Two sample Kolmogorov-Smirnov test
						</a>
					</li>
					
					<li>
						<a href="/20190220/threejs-test-page.html" title="three.js test page">
								<img class="small-icon" src="/assets/img/tree.png" title="three.js test page" />three.js test page
						</a>
					</li>
					
					<li>
						<a href="/20181111/hic_data_format.html" title="Hi-C 数据分析结果应该怎么存？">
								<img class="small-icon" src="/assets/img/file.ico" title="Hi-C 数据分析结果应该怎么存？" />Hi-C 数据分析结果应该怎么存？
						</a>
					</li>
					
					<li>
						<a href="/20181010/d3_bubble_chart.html" title="用 D3.js 画一个 bubble chart">
								<img class="small-icon" src="/assets/img/file.ico" title="用 D3.js 画一个 bubble chart" />用 D3.js 画一个 bubble chart
						</a>
					</li>
					
					<li>
						<a href="/20180724/new_kind_slides.html" title="论制作 Slides 的几种姿势">
								<img class="small-icon" src="/assets/img/file.ico" title="论制作 Slides 的几种姿势" />论制作 Slides 的几种姿势
						</a>
					</li>
					
					<li>
						<a href="/slides/test/slideshow.html" title="Markdown Slides Test">
								<img class="small-icon" src="/assets/img/slides.png" title="Markdown Slides Test" />Markdown Slides Test
						</a>
					</li>
					
					<li>
						<a href="/20180221/bioview.html" title="bioView - 一个生信常用文件格式的可读性增强工具">
								<img class="small-icon" src="/assets/img/file.ico" title="bioView - 一个生信常用文件格式的可读性增强工具" />bioView - 一个生信常用文件格式的可读性增强工具
						</a>
					</li>
					
					<li>
						<a href="/20170916/markdown-test-page.html" title="Markdown Test Page">
								<img class="small-icon" src="/assets/img/file.ico" title="Markdown Test Page" />Markdown Test Page
						</a>
					</li>
					
					<li>
						<a href="/20170831/gol-js.html" title="一个JS实现的生命游戏">
								<img class="small-icon" src="/assets/img/file.ico" title="一个JS实现的生命游戏" />一个JS实现的生命游戏
						</a>
					</li>
					
					<li>
						<a href="/20170831/parallel.html" title="这大概是程序串行改并行最简单粗暴的方法">
								<img class="small-icon" src="/assets/img/file.ico" title="这大概是程序串行改并行最简单粗暴的方法" />这大概是程序串行改并行最简单粗暴的方法
						</a>
					</li>
					
					<li>
						<a href="/20170730/learn-docker.html" title="学习Docker">
								<img class="small-icon" src="/assets/img/file.ico" title="学习Docker" />学习Docker
						</a>
					</li>
					
					<li>
						<a href="/20170530/learn-assemble.html" title="学习汇编语言">
								<img class="small-icon" src="/assets/img/file.ico" title="学习汇编语言" />学习汇编语言
						</a>
					</li>
					
					<li>
						<a href="/20170505/schoolnet.html" title="如何在非校园网环境下使用学校文献数据库">
								<img class="small-icon" src="/assets/img/file.ico" title="如何在非校园网环境下使用学校文献数据库" />如何在非校园网环境下使用学校文献数据库
						</a>
					</li>
					
					<li>
						<a href="/20170308/zerotier.html" title="使用ZeroTier搭建虚拟局域网">
								<img class="small-icon" src="/assets/img/file.ico" title="使用ZeroTier搭建虚拟局域网" />使用ZeroTier搭建虚拟局域网
						</a>
					</li>
					
					<li>
						<a href="/20161224/hy-in-brief.html" title="Python生态下的Lisp方言">
								<img class="small-icon" src="/assets/img/file.ico" title="Python生态下的Lisp方言" />Python生态下的Lisp方言
						</a>
					</li>
					
					<li>
						<a href="/20161210/scrapy_douban.html" title="使用Scrapy爬取豆瓣相册">
								<img class="small-icon" src="/assets/img/file.ico" title="使用Scrapy爬取豆瓣相册" />使用Scrapy爬取豆瓣相册
						</a>
					</li>
					
					<li>
						<a href="/20161115/speed_up_python.html" title="加速Python">
								<img class="small-icon" src="/assets/img/file.ico" title="加速Python" />加速Python
						</a>
					</li>
					
				</ul>
			
		</div>
		<div class="post_total">
			
				<div class="left">25 object(s)</div>
			
			<div class="right">&nbsp;</div>
		</div>
	</div>
	
        <div class="content">
			<div class="post_title">
				<img src="/assets/img/file.png" />
				<h1>加速Python</h1>
				<a href="/"><div class="btn"><span class="fa fa-times"></span></div></a>
				<div class="btn btn_max"><span class="fa fa-window-maximize"></span></div>
				<div class="btn"><span class="fa fa-window-minimize"></span></div>
			</div>
			<ul class="topbar">
				<li>November 15, 2016</li>
			</ul>
			<div class="post_content" style="max-height: 600px">
				<div class="post_content_inner">
        		<p>C++的作者B.Stroustrup 说过世界上只有两种语言，一种经常有人对它抱怨，还有一种根本没人用(比如___，免得我被打请大家自行填空_(┐「ε:)_
)。</p>

<p>Python最经常被人抱怨的就是慢（大概比c要慢40倍，的确是挺感人的…具体数据可以参考<a href="http://benchmarksgame.alioth.debian.org/u64q/compare.php?lang=python3&amp;lang2=gcc">这里</a>）和它的<a href="https://wiki.python.org/moin/GlobalInterpreterLock">gil</a>（没有真正的线程这一点是经常遭到诟病的地方）。的确，现在cpu已经很快很快了，对于Python这样一门表达能力超强的语言，我们在一定程度上可以对它的速度有所忍受，并且对于网络应用来说python通过异步来实现并发也是可以非常快的。但Python 有时也会被用来做一些计算密集型的工作，这个时候就需要一些办法来加速了，现在看来有一大票办法来加速Python，并且其中有一些能够一定程度上解决gil的问题，最近我学习了其中几种，这里就结合一个例子来简单谈一下感受。</p>

<!-- more -->

<h2 id="一个计算密集型过程例子julia-set">一个计算密集型过程例子：<a href="https://en.wikipedia.org/wiki/Julia_set">julia set</a></h2>

<p>julia set 与 <a href="https://en.wikipedia.org/wiki/Mandelbrot_set">mandelbrot set</a> 一样，是一种<a href="https://en.wikipedia.org/wiki/Fractal">分形结构</a>，它能由极其简单的规则所描述，但产生出的图形却十分的迷人与复杂。julia set的计算规则真的非常简单，以便理解后面的代码我先说明一下，已经理解的朋友请跳过这一部分。
julia set可以由<img src="https://wikimedia.org/api/rest_v1/media/math/render/svg/191627a3eebdd6608c9b226786defc468b747502" alt="" />这个复数迭代公式来描述，其中z和c都是复数，现在我们把复数空间对应到一个二维平面上，比如实数部分对应于x轴，虚数部分对应y轴，然后如果我们固定一个c的值并且对这个二维中的每一个点使用上面的公式进行迭代，会发现有的点经多次迭代后会收敛到某一范围内，另一些则会无限增大，如果我们把所有收敛的点在这个二维平面上画出来，我们将会的到非常奇妙的分形图案，并且改变c的值，图案也会发生变化，相关的图片可以看看julia set的<a href="https://en.wikipedia.org/wiki/Julia_set">维基</a>页面。</p>

<p><img src="https://upload.wikimedia.org/wikipedia/commons/b/b1/Julia_set_%28ice%29.png" alt="julia set" width="95%" /></p>

<p>现在我们来尝试用python来计算出一个julia set的图形，并用几种方法对这个过程进行加速，来看看它们的加速效果。</p>

<h3 id="pure-python">pure Python</h3>

<p>首先我们用纯Python实现julia set的计算，为对比后面各种手段的加速效果提供参考。</p>
<pre><code class="language-python"># jl.py 
# calc_julia pure python edition

import numpy as np

def escape(z, c, z_max, n_max):
    i = 0
    while (z.real**2 + z.imag**2 &lt; z_max**2) and i &lt; n_max:
        z = z*z + c
        i += 1
    return i

def calc_julia(resolution, c, bound=1.5, z_max=4, n_max=1000):
    step = 2.0 * bound / resolution
    counts = np.zeros((resolution + 1, resolution + 1), dtype=np.int32)
    for i in range(resolution + 1):
        real = -bound + i * step
        for j in range(resolution + 1):
            imag = -bound + j * step
            z = real + imag * 1j
            counts[i, j] = escape(z, c, z_max, n_max)
    return np.asarray(counts)
</code></pre>
<p>注意这里虽然import了numpy但是并没有使用numpy进行计算，只是把计算结果存储在一个numpy数组里方便以使用matplotlib进行绘图。然后我们可以使用calc_julia这个函数计算得到julia set并对它进行可视化了。
打开ipython notebook：</p>
<pre><code class="language-python">In [1]: import matplotlib.pyplot as plt
In [2]: import numpy as np
In [3]: from jl import calc_julia
In [4]: jls = calc_julia(1000, -0.4+0.6j)
In [5]: plt.imshow(np.log(jls));plt.show()
</code></pre>
<p>这样我们就会得到对于复数 -0.4 + 0.6i 的julia set 图形，如下，感觉挺神奇的╰( ᐖ╰)。</p>

<p><img src="/images/blog/julia_set_-0.40.6i.png" alt="julia_set_-0-40-6i" /></p>

<p>这很有趣但是计算jls这一步的确是慢了一点，大概要7.6秒，如果我们想多看几个其他复数的julia set，恐怕是等的会有点烦了。我们下面来看一看加速后这个过程要多久吧。</p>

<h3 id="numba">numba</h3>

<p><a href="http://numba.pydata.org/">numba</a>大概是使用起来最简单的计算加速办法了吧，只需要在函数上面加一个小小的装饰器，numba就会把这个函数编译成<a href="http://llvm.org/">LLVM</a>的IR码从而达到<a href="https://en.wikipedia.org/wiki/Just-in-time_compilation">jit</a>的加速效果。我们来试一下。</p>
<pre><code class="language-python"># jl.py
# calc_julia numba edition

import numpy as np
from numba import jit

@jit
def escape(z, c, z_max, n_max):
    i = 0
    while (z.real**2 + z.imag**2 &lt; z_max**2) and i &lt; n_max:
        z = z*z + c
        i += 1
    return i

@jit
def calc_julia(resolution, c, bound=1.5, z_max=4, n_max=1000):
    step = 2.0 * bound / resolution
    counts = np.zeros((resolution + 1, resolution + 1), dtype=np.int32)
    for i in range(resolution + 1):
        real = -bound + i * step
        for j in range(resolution + 1):
            imag = -bound + j * step
            z = real + imag * 1j
            counts[i, j] = escape(z, c, z_max, n_max)
    return np.asarray(counts)
</code></pre>
<p>看很简单吧，只需要在函数上面加上@jit这个装饰器即可。测试一下速度，完成与上面相同的计算只需要0.6s比纯python快了12.6倍左右。</p>

<h3 id="pypy">PyPy</h3>

<p><a href="http://pypy.org/features.html">pypy</a>是一个jit的Python解释器实现，不需要修改代码就能执行，也相对比较方便，但是毕竟它已经脱离了CPython的生态，配置起来比起numba那样直接一条pip指令就能安装还是稍微麻烦一点，并且还有好多Cython中的库在这里也是没有的。把上面纯Python的代码用PyPy运行一遍,耗时0.44s，效果很可观，甚至比numba还要快一点。</p>

<h3 id="cython">Cython</h3>

<p>Cython即带有C类型的Python，通过引入C类型把Cython代码编译成高效的C代码，相当于使用类型信息对Python代码的“蒸馏”，并且“蒸馏”的粒度是可控的，纯Python代码在Cython中完全是合法的，只是引入越多的类型信息编译出来的代码就越高效，并且还能用来方便地包装C代码和使用C库，可以说Cython是这些方案之中最强大最成熟的一种，很多大型项目诸如Scipy、Pandas、sklearn中都大量使用了Cython。它用来写库是一种非常好的选择，对于Python用户来说Cython是非常值得学习的，如果感兴趣的话，有一本很好的讲Cython的书<a href="http://shop.oreilly.com/product/0636920033431.do">Cython A Guide for Python Programmers</a>可以找来阅读一下，网上很容易找到电子版。以下的计算julia set的Cython代码就来源于这本书。</p>
<pre><code class="language-python"># jl.pyx
# calc_julia cython edition

from cython import boundscheck, wraparound
import numpy as np

cdef int escape(double complex z,
                double complex c,
                double z_max,
                int n_max) nogil:
    cdef:
        int i = 0
        double z_max2 = z_max * z_max
    while norm2(z) &lt; z_max2 and i &lt; n_max:
        z = z * z + c
        i += 1
    return i

cdef inline double norm2(double complex z) nogil:
    return z.real * z.real + z.imag * z.imag

@boundscheck(False)
@wraparound(False)
def calc_julia(int resolution, double complex c,
                double bound=1.5, double z_max=4.0, int n_max=1000):
    cdef:
        double step = 2.0 * bound / resolution
        int i, j
        double complex z
        double real, imag
        int[:, ::1] counts
    counts = np.zeros((resolution + 1, resolution + 1), dtype=np.int32)
    for i in range(resolution + 1):
        real = -bound + i * step
        for j in range(resolution + 1):
            imag = -bound + j * step
            z = real + imag * 1j
            counts[i, j] = escape(z, c, z_max, n_max)
    return np.asarray(counts)
</code></pre>
<p>恩比起numba与pypy似乎麻烦了一点，多了很多以cdef开头的类型声明，不过当你掌握了它简单的语法以后会发现这并不难就像Python一样。编写完.pyx文件后与之前不同的是这里需要编译Cython代码，所以我们再创建一个名为setup.py的文件。</p>
<pre><code class="language-python"># setup.py
from distutils.core import setup
from Cython.Build import cythonize

setup(name="julia",
    ext_modules=cythonize("jl.pyx"))
</code></pre>
<p>然后在终端内敲入python setup.py build_ext -i进行编译，然后就能与之前一样的方式使用编译好的函数了，运行一下，计算过程耗时0.16s，比pypy更快了，比最初的纯Python版本快了40多倍！</p>

<h2 id="结语">结语</h2>

<p>这里说了三种加速Python的办法，总结一下就是写独立的应用(不需要太多的第三方库的)可以试试PyPy，加速函数运行速度就用numba，写库就用Cython。其实还有很多加速方案啦，比如使用将计算向量化使用numpy,以及它的多线程版numexper，以及使用多进程multiprocessing多进程编程等等，具体选哪种还要看情况，当然对于优化速度来说最重要的还是选对合适的算法。并且在优化之前要想清楚现阶段需不需要优化，很多时候对于比如科学计算这样的应用，比起速度可能更在意的还是正确性，先把功能做对了再来看优化吧。</p>

				
					<br>
<hr>
<br>
<div class="donate">
	<p>Email: nanguage@yahoo.com</p>
</div>
				
				</div>
			</div>
		</div>
    
	<script src="/assets/js/001.js"></script>
	<script src="/assets/js/002.js"></script>
	<div class="footer">
		<p>blog theme: <a href="https://github.com/h01000110/windows-95">win95</a></p>
	</div>
</body>
</html>